Magnons and magnetic fluctuations in atomically thin MnBi2Te4

Electron band topology is combined with intrinsic magnetic orders in MnBi2Te4, leading to novel quantum phases. Here we investigate collective spin excitations (i.e. magnons) and spin fluctuations in atomically thin MnBi2Te4 flakes using Raman spectroscopy. In a two-septuple layer with non-trivial topology, magnon characteristics evolve as an external magnetic field tunes the ground state through three ordered phases: antiferromagnet, canted antiferromagnet, and ferromagnet. The Raman selection rules are determined by both the crystal symmetry and magnetic order while the magnon energy is determined by different interaction terms. Using non-interacting spin-wave theory, we extract the spin-wave gap at zero magnetic field, an anisotropy energy, and interlayer exchange in bilayers. We also find magnetic fluctuations increase with reduced thickness, which may contribute to a less robust magnetic order in single layers.


H-field dependent Raman spectra
Systemic field dependent Raman spectra from -20 cm −1 to 20 cm −1 without subtracting the spin fluctuations background discuss in the main text. The left(right) panel shows spectra observed in the co-circular(cross-circular) polarization channel. Figure S2: Magnetic field dependent Raman spectra from 2-SL MBT. (a) Low temperature (12 K) co-circular (σ + /σ + ) polarization spectra. The blue, green, and red lines highlight the magnon peak in the three magnetic phases discussed in the main text. (b) Low temperature (12 K) cross-circular (σ + /σ − ) polarization spectra. In both (a) and (b), the shaded gray rectangles block a noise peak that originates from the instrument.

Magnetic lattice and Effective spin model
The magnetic lattice structure for 2-SL MBT is shown in Fig. S3. In each layer, the Mn atoms form a triangular lattice in a hexagonal setting with axes a = A(1/2, − √ 3/2, 0), b = A(1/2,  The six (three) intralayer (interlayer) nearest-neighbors are shown by black (red) lines.
The local-moment Hamiltonian (S = 5/2) deduced from inelastic neutron scattering mea-surements 1 for bulk MBT in the presence of an applied external magnetic field is where J ij describes intralayer interactions (up to fourth-order neighbor interactions are need to correctly fit the neutron scattering data with SJ 1 = 0.3 meV, SJ 2 = −0.083 meV, and SJ 4 = 0.023 meV), SJ c = −0.055 meV corresponds nearest-neighbor AFM interlayer interactions, and In this work, we adopt the model Hamiltonian Eq. (1), but restrict the intralayer interactions to nearest-neighbors. We use the interlayer exchange interaction J c and single-ion anisotropy D as fitting parameters to describe our Raman scattering measurements.

Spin wave calculations
In this section, we calculate the spin wave spectrum for 2-SL in the AFM, canted AFM (cAFM), and FM phases as a function of magnetic field. Since the primitive unit cell contains two magnetic atoms in all the magnetic phases, we expect two magnon modes across all the magnetic fields considered.
Spin order as a function of magnetic field. 2-SL is expected to exhibit interlayer AFM order in the absence of applied magnetic field. When a magnetic field is applied in the c-direction, the system undergoes a first-order spin-flop phase transition into a cAFM order at the critical field H sf 3 . Upon further increasing the field, the system becomes FM at a second critical field H f m . In Table 1, we summarize the reported values for the cAFM and FM critical fields for several samples.
To determine the ground state in our sample for a given magnetic field, we consider the macro-spin approximation 5 , where the energy per unit volume derived from the local moment Hamiltonian (1) takes the form where we assume M i = M s (sin θ i , 0, cos θ i ), M s is the saturation magnetization per unit cell, z 1 is the number of intralayer n.n.s, and z c is the number of interlayer n.n.s. in the 2-SL system. Solutions to the set of equations ∂E(M 1 , M 2 )/∂θ i = 0 determine the possible ground states. The AFM is given by θ 1 = 0 and θ 2 = π, while the spin-flop phase by θ 1 = −θ 2 = arccos (−gµ B H z /(2z c SJ c + 2DS)).
AFM magnons When the system is in an AFM groundstate, we perform the Holstein-Primakoff . We keep terms in the transformation up to quadratic order in the bosonic operators. After Fourier transforming to momentum space we obtain the Bogoliubov-de Gennes (BdG) where and Diagonalizing H(k) we obtain the magnon energies In Fig. S4, we plot the Γ point AFM magnon energies as a function of the applied magnetic field. The corresponding eigenvectors at the Γ point are where N is a normalization constant. These wavefunctions indicate that the two sublattices oscillate with different amplitudes. This has consequence for the selection rules, as discussed in the main text.
In Fig. S5, we plot the magnon energies for H z = 0 T along a high-symmetry path in the BZ. The red bands correspond to the phonons in the AFM phase obtained with density functional theory as outlined in Note 6. Notice that in the vicinity of the Γ point, the mangnons overlap with the acoustic phonons.
FM magnons For the FM case, which applies for magnetic fields above the saturation field H F M , we apply a Holstein-Primakoff transformation to the local moment Hamiltonian (1), defined bŷ H sf where k = (k x , k y ) is the momentum defined in the Brillouin zone (BZ) of the triangular lattice.
The magnon energies are At the Γ point, ω ± (Γ) = 2SD + gµ B H z + z c S(J c ± J c ). The intralayer exchange interaction does not enter in the energy at the Γ point. The associated wave functions are Therefore, the ω + (Γ) branch corresponds to magnetic moments in opposite layers oscillating in opposite directions, while the ω − (Γ) branch corresponds to sublattices oscillating in the same direction.
In Fig. S6, we plot the magnon energies as a function of the applied magnetic field. In Fig. S7, we plot the magnon energies along a high-symmetry path in the BZ. We also display the phonon energies (red curves).
Canted AFM magnons In this section we consider the magnons in the canted AFM phase (cAFM), obtained in the magnetic field regime H sf < H z < H f m T. The canting angle is given by . For this derivation, we follow closely Ref. 5 . We introduce a coordinate system for each layer, where the local z-axis is aligned with the canted moments. We apply rotations R y (θ) and R y (−θ) to each sublattice. Then, we apply a Holstein-Primakoff transformation to the rotated spin operatorsŜ where  At the Γ point, the cAFM magnon energies simplify to In Fig. S8, we plot the cAFM magnon energies as a function of the magnetic field.

Magnon selection rules
In this section, we discuss the selection rules for the magnon modes in 2-SL MBT in each of the three ground states as a function of magnetic field. In analogy with the selection rule derivation for phonons, we need to take into account the magnetic point group and analyze the symmetries of the magnetic moment oscillations. For our symmetry analysis, we employ the ISOTROPY Software Suite 6, 7 , and the Bilbao Crystallographic Server 8 .
In the paramagnetic phase, 2-SL MBT belongs to the space group P3m1 with point group We can understand this difference in space group compared with bulk MBT because 2-SL MBT (and N-SL in general) loses the symmetries composed with translations along the c axis. In 2-SL MBT, neither FM nor AFM order changes the size of the chemical unit cell, and the wave vector for the magnetic order is k = (0, 0, 0). This is relevant to derive the possible symmetry-allowed magnetic orders. For an extended discussion, see Ref. 9 . In the next sections, we obtain the magnetic point group in each phase, the corresponding Raman tensors, and the expected one-magnon selection rules. We then study two-magnon Raman processes, and the effect of broken inversion symmetry.
The corresponding magnetic point group is where T is the time reversal operator, i corresponds to inversion, and the point group  Table 2.
According to the group theory of magnetic point groups, the selection rules for the Γ point where all the coefficients are real. The Raman intensity is then given by I ν ∝ α,β={x,y,z} e * ponents of the Raman tensor for mode ν 15 . Then, a magnon with symmetry DA 1 can be observed in the σ + /σ + channel, and magnons with symmetry DE can be observed in the σ + /σ − channel.
The corepresentation of the magnons in AFM 2SL-MBT can be obtained by studying the transformation properties of the magnon wave functions derived in the previous sections 16 . At the Γ point, we find |ψ − (k = 0) = 1/N Câ +b † and |ψ + (k = 0) = 1/N Cb +â † , where C is a coefficient that depends on the interlayer exchange and anisotropy (see Eqns. (7)).
cAFM phase In the cAFM phase, the magnetic space group is C2 , as determined by ISOTROPY 6, 7 .
The corresponding not centro-symmetric magnetic point group 2 = {1, T 2 y }. The Raman tensor for the magnetic point group 2 corepresentations is 14 In this magnetic phase, due to the low-symmetry, the Raman tensor components do not have many restrictions, and a magnon could be observed in the σ + /σ + and σ + /σ − polarization channels, depending on the intensity of the matrix elements. These intensities cannot be obtained within group theory and require first principles calculations.
Spin wave Raman tensors We list the Raman elements of all possible observable modes in different magnetic phases in MBT 2-SLs. Experimentally, not all modes allowed by the Raman tensor are observed. One needs to compute the magnon wavefunctions for a particular material using spin wave theory. The calculated magnon modes obey a certain symmetry (e.g. DA 1 , DE) that further predicts the polarization configurations of the Raman signal.

Magnetic
Raman

Symmetry Operations
The c-AFM phase belongs to the 2' point group where two-fold rotational symmetry along the y-axis ( C y ) is broken. On the other hand, in the AFM case two-fold rotational symmetry along the xy plane (C y ,C yx ,C x ) are lost. So in the c-AFM case, one must rotate along the y-axis while the AFM case can be either x or y-axis. Figure S9: (a -c) Illustration of the step by step symmetry operations for each magnetic phase of 2-SL MBT.
In summary, even though the symmetries are different in each of the magnetic ground states due to the direction of the magnetic moments, one-magnon processes can be observed in the σ + /σ + polarization channel for cAFM and FM phases. In the AFM phase, the one-magnon processes are expected in the cross-polarization channel σ + /σ − .

Two-magnon Raman scattering
In this section, we consider the possibility that the magnetic mode at zero applied magnetic field originates from two-magnon Raman scattering processes.
There are two features that support this interpretation. To model two-magnon Raman scattering processes, we consider the Loudon-Fleury scattering operator for spin-zero excitations [21][22][23][24] to opposite sublattices (in our case, opposite layers), andê in(out) is the incoming (outgoing) electric field polarization. In the σ + /σ + channel, we find that the Raman intensity 22,25,26 , where |0 is the ground state with energy ω 0 , {|µ , ω µ } are the two-magnon energy and eigenstates obtained from the spin-wave Hamiltonian in the noninteracting limit. We find that the two-magnon Raman intensity I(ω) captures all the characteristics of the density of states, as shown in Figure S11, where we consider the interlayer exchange interaction in the Loudon-Fleury scattering operator. Since this is a spin zero excitation, it has no magnetic field dependence.
Therefore, the absence of magnetic field dependence of a zero-spin excitation and the energy scale support the hypothesis that the magnetic mode at zero and low-applied magnetic field can be associated to a two-magnon process.
We note that in LiMnPO 4 , an antiferromagnetic material with a spin-flop transition at ∼ 4 T 26 , the two-magnon responses are observed for magnetic fields above the transition field 25 .   Figure S11: Two magnon Raman intensity as a function of energy in the σ + /σ + channel. The red arrow indicates the step-like feature possibly associated with experimental observation.

Effect of broken inversion symmetry on the Raman selection rules
In the paramagnetic phase, the point group of the crystal is D 3d , which possess inversion symmetry. When inversion symmetry is broken, for example as the result of strong coupling with the substrate, the point group reduces to C 3v . The correlation relation table is 8 We note that below the transition temperature AFM order sets in, which also breaks inversion symmetry, and can also change the phonon selection rules. The phonon spectrum is presented in Fig. S14. In Table 6, we show the phonons at the Γ point, along with the irreps.   Figure S14: Phonon spectrum along a high-symmetry path in the Brillouin zone.

Magnetic fluctuations in paramagnetic phase
Quasielastic scattering In this section, we discuss quasielastic scattering (QES) in paramagnetic phase. We first apply Bose factor correction to obtain normalized Raman susceptibility: where n is the Bose-Einstein factor and I(ω) is the Raman intensity. We observe clear background changes as a function of the layer thickness (Fig. S15). We attribute this background to quasielastic scattering for magnetic fluctuations. The spin fluctuations couple the light to the magnetic energy, giving rise to the peak at zero Raman frequency shift, the linewidth of spin-lattice relaxation time, and the intensity that is proportional to magnetic specific heat. We extract the quasielastic scattering using a magnetic fluctuation model, where C m is magnetic specific heat and D is the magnetic contribution to thermal conductivity [33][34][35] . The quasielastic scattering is enhanced in thinner layers. This illustrates that the magnetic fluctuations are stabilized as the magnetic system approaches to bulk regime. Moreover, the spin fluctuations involve the relaxation channel by coupling with lattice vibrations 33,36 . The spin-lattice coupling was observed in our previous study.
8 Temperature dependent dynamic Raman susceptibility Figure S16: Temperature dependence of dynamic Raman susceptibility for 2 and 3 SL. The grey shade indicates the critical behavior across the Nèel temperature.
We analyze temperature dependence of dynamic Raman susceptibility χ dyn in 2 and 3 SLs (Supplementary Figure S17). We observe χ dyn is peaked at the Nèel temperature, ∼ 17 K. The observed Nèel temperature for 2 SLs is consistent with the integrated intensity analysis shown in Fig. 2c. The difference between 2 and 3 SLs is unclear due to the limited temperature step (2 K).
χ dyn is proportional to magnetic specific heat C m , which shows the critical behavior of magnetic phase transition 34, 35 . This temperature dependence support our analysis of the QES as evidence for increasing magnetic fluctuations with decreasing layer thickness.

SL 2 SL 3 SL 4 SL 5 SL
Intensity (a.u.) Raman Shift (cm -1 ) Figure S17: Layer dependence of low frequency Raman spectra in 1 to 5 SLs. The data is taken at room temperature. The stripe-pattered shade is to block the noise line. The solid lines are Lorentzian model fits.
Low-frequency phonon spectra taken at room temperature using collinearly polarized incident and scattered light. The phonon mode corresponds to a breathing mode. The systematic frequency shift as a function of layer thickness can be modeled by a linear chain model. The mode observed from the 1-SL corresponds to lattice vibration against the substrate. The presence of this mode suggest the crystalline structure of 1-SL is retained.